function [allsols,  scrs, bestslns] = nonlinearOptimizeselfcalibnormMODTRACE(TF,w,h,numParams,initialWeights,numTries,kstart)

numF=size(TF,2);

if (nargin <5)
    initialWeights=ones(numF,1);
end
if (nargin <6)
    numTries=1;
end
if (nargin <7)
    kstart=[w 0 w/2; 0 w h/2; 0 0 1];
end



allsols= repmat({zeros(3)}, 1, numTries);
scrs=zeros(numTries,1);
bestslns=zeros(3,3);
initialSolutions=createInitialSolutions(numTries,w,h,numParams,kstart);
[lb,ub]=findBoundsonK(w,h,numParams);
bestscore=100000000;


for i=1:numTries
    
    [x,error]=nonlincustomFunction(initialSolutions{1,i},TF,initialWeights,w,h,lb,ub);
    
    allsols{1,i}=convertXTOKselfK(x,w,h);
    % display([ num2str(i) ' out of ' num2str(numTries)]);
    scrs(i,1)=error;
    
    if(scrs(i,1)<bestscore  )
        
        bestscore=scrs(i,1);
        bestslns=allsols{1,i};
        
    end
    
end
end
function [x,error]=nonlincustomFunction(x0,F,weights,w,h,lb,ub)
% according to the discussion weights need to be cholesky decomposed
CholeskyinitialWeights = sqrt(weights);
CholeskyinitialWeights=1000*CholeskyinitialWeights*(1/sum(CholeskyinitialWeights));

levmaroptions=[1E-5, 1E-25, 1E-25, 1E-30, 1E-06]; % this is from the demos


sizeFs=size(F,2);
numParams=size(x0,2);

if(numParams>sizeFs)
    measurements= zeros(numParams,1);
    outputformat=2 ;
else
    outputformat=1  ;
    measurements= zeros(sizeFs,1);
end
type=1;
if(type==1)
    [ret, x, info]=levmar('computerEssentialErrorSVDNFramesGeneralMODTRACE','computerEssentialJACOBIANSVDTRACE' ,  x0, measurements, 2000, levmaroptions,'bc',lb, ub, F,CholeskyinitialWeights,w,h,outputformat);
    
    error=sum(abs(info(2)));
else
    outputformat=3;
    tolx=1e-41; tolf=1e-41;
    optionsfmincon = optimset('GradObj','on', 'Algorithm','interior-point','Display','on','TolFun',tolf,'TolX',tolx);
    f = @(x)computerEssentialErrorSVDNFramesGeneralMODTRACE(x,F,CholeskyinitialWeights,w,h,outputformat);
    
    %'DerivativeCheck', 'on',
    
    try
        [x,fval,exitflag,output] = fmincon(f,x0,[],[],[],[],lb,ub,[],optionsfmincon);
    catch err
        display('*********** caught some error');
        x=x0;
        fval=10;
        save('badcrashfmincon.mat','f','x0','lb','ub','optionsfmincon');
    end
    error=sum(fval);
end

end

function [initialSolutions]=createInitialSolutions(numTries,w,h,numparams,KSTART)

K=KSTART;

fstd=400;
skewmax=12;
initialSolutions=repmat({zeros(1,numparams)}, 1, numTries);


for i=1:numTries
    Knew=K;
    if(i>1)
        Knew(1,1)=initialSolutions{1,1}(1,1)+randn()*fstd;% gaussian
        Knew(2,2)=initialSolutions{1,1}(1,1)+randn()*fstd;% gaussian
        Knew(1,3)=rand()*w;
        Knew(2,3)=rand()*h;
        Knew(1,2)=rand()*skewmax;
        
        
    end
    initialSolutions{1,i}=convertKTOXselfK(abs(Knew),numparams);
end

end

function [lb,ub]=findBoundsonK(w,h,numParams)
fmax=2000;  fmin=40;
skewmax=10; skewmin=-1;
xmax=w;     xmin=0;
ymax=h;     ymin=0;
armax=1.2;  armin=0.8;
KMAX = [fmax   skewmax        xmax; 0      fmax*armax   ymax;  0            0             1  ];
KMIN = [fmin   skewmin        xmin; 0      fmin*armin   ymin;  0            0             1  ];

lb= convertKTOXselfK(KMIN ,numParams);
ub= convertKTOXselfK( KMAX ,numParams);

end